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ABSTRACT 

We calculate the expected nHz-/iHz gravitational wave (GW) spectrum from coalescing Mas- 
sive Black Hole (MBH) binaries resulting from mergers of their host galaxies. We consider detec- 
tion of this spectrum by precision pulsar timing and a future Pulsar Timing Array. The spectrum 
depends on the merger rate of massive galaxies, the demographics of MBHs at low and high red- 
shift, and the dynamics of MBH binaries. We apply recent theoretical and observational work 
on all of these fronts. The spectrum has a characteristic strain hc{f) ~ 10~^^(//yr~^)~^/^, just 
below the detection limit from recent analysis of precision pulsar timing measurements. However, 
the amplitude of the spectrum is still very uncertain owing to approximations in the theoretical 
formulation of the model, to our lack of knowledge of the merger rate and MBH population at 
high redshift, and to the dynamical problem of removing enough angular momentum from the 
MBH binary to reach a GW-dominated regime. 

Subject headings: 

1. Introduction 

Over the past decade, evidence has accumulated that Massive Black Holes^ (MBHs) are ubiquitous 
in the spheroidal components of low-redshift galaxies {e.g., Magorrian et al. 1998), and that their masses 
are tightly correlated with the properties of the host galaxies. These initial studies noted that the MBH 
mass scaled with the mass of the spheroidal component of the parent galaxy as inferred from the luminosity 
although there was considerable scatter. More recently, a much tighter correlation has been observed between 
the MBH mass and the velocity dispersion in the spheroid (Ferrarese & Merritt 2000; Gebhardt et al. 2000). 

At higher redshift, evidence for the ubiquity of MBHs is more circumstantial. However, we do know 
that many distant galaxies harbor active nuclei for some portion of their life and, in turn, that the engine 
driving this activity is almost certainly gravitational accretion onto a central black hole. There are many 
outstanding questions regarding the era and duty cycle of AGN activity, the initial mass and mass evolution 
of MBHs, and the mass growth of MBHs during AGN episodes. Nonetheless, the existing evidence strongly 
favors the hypothesis that the present-day MBH population is the dormant remnant of the AGN population 
at high redshift. 



^The literature has variable usage of the terms "Massive" and "Supermassive" to describe the 10® ^ Mq black holes in the 
centers of galaxies (e.g., Hug ????). For concision, we favor the former in this work. 
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During the last decade astrophysicists have also begun to understand the important role that mergers 
play in the history of galaxies. In the hierarchical clustering picture, mergers are by definition responsible 
for building up the mass of today's galaxies. Numerical simulations have also shown that the merger history 
also determines the morphology of present-day galaxies: e.g., major mergers tend to produce elliptical 
galaxies. The strongest evidence for the mergers is the relatively flat density scaling laws in "some" spheroids 
(Milosavljevic & Merritt 2001). Simulations that reproduce this scaling actually calculate the merger of 
galaxy halos where the majority of the mass resides. The subsequent merger of the central galaxies within 
the halos is nearly guaranteed owing to the rapid process of dynamical friction for these dynamically soft 
stellar systems. 

The ubiquity of MBHs in galaxies and the frequent growth of galaxies by merger leads to the obvious 
question about the fate of the pair of MBHs in the subsequent evolution of merged systems. The same dy- 
namical friction process that brings together parent galax;ies within merged halos will also lead to the sinking 
of the MBHs to the center of the merger daughter galaxy. Unfortunately, in the simplest case dynamical 
friction will "turn off" well before the MBH binary coalesces^, or even before it reaches a regime where 
gravitational radiation losses will drive the binary's evolution (Begelman et al. 1980). Various mechanisms 
have been suggested to evolve the MBH pair through this period. 

If the pair can indeed evolve to the gravitational wave (GW) regime, it will over time lose angular 
momentum to gravitational radiation and eventually coalesce. The final inspirals of MBH pairs take about 

10** seconds, and arc the most luminous gravitational wave events in the Universe. These events are one 
of the chief targets of the Laser Interferometer Space Antenna (LISA) GW interferometer satellite. In this 
paper, we concentrate on the quiescent evolution leading up to this state. Before the final inspiral, the GW 
amplitude is 

4.4 X IQ-''MI''P--I^D-^1^^^, (1) 

where h is the strain or metric perturbation, Mg is the total mass of the binary in units of 10* M©, q is the 

mass ratio {q < 1), Pyr is the observed GW period in years (the orbital period divided by 2 for the expected 
nearly-circular orbits), and i?Gpc is the distance in Gpc. The lifetime of the system is 

tGW = 1.1 X 10^ yr M8-^/^P«/3ii±^. (2) 

q 

The precision of the rotation pcaiods of millisecond period pulsars which is established via pulse arrival 
time measurements allows detection of the stochastic MBH-MBH background spectrum at nHz frequencies 
(Sazhin 1978; Detweiler 1979; Rajagopal & Romani 1995). The detection process is very similar to that of 
the laser interferometers — passage of an electromagnetic signal through distorted space-time — except that 
in the case of pulsars no mirrors are needed as we have a distant precision clock, the rotating neutron star, 
sending us a periodic pulse train. Measurements of a single pulsar can place an upper limit on the presence 
of a spectrum of gravitational radiation (Kaspi et al. 1994; Lommen 2001). In short, the limit on strain is 
given by the ratio of the timing precision and the measurement duration; e.g. 1 /is/10 yr ~ 3 x 10~^^. 
While pulsar timing cannot detect individual events of the small amplitude given in equation (1), these 
measurements have a good prospect for detection of the stochastic background of GWs from the ensemble 



^We use the term merger to denote the joining of halo and galaxy pairs and the term coalescence for the joining of the 
members of a MBH binary to distinguish implicitly the two nominally sequential processes. 
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of these MBH-MBH coalescence events throughout the Universe^. The GW perturbs pulse arrival times of 
a spatial array of pulsars in a correlated manner that is distinct from other known perturbations such as 
atomic time and ephemeris errors. A Pulsar Timing Array then acts as a nHz gravitational wave telescope 
capable of direct detection of the stochastic background radiation. 

Questions remain. Are there dynamical processes that drive the MBH pair into the GW regime? Does 
this happen frequently? If so, can we observe the gravitational radiation via pulsar timing? If not, do 
we then see evidence of close binary MBHs "hung up" in the centers of many, or most, galaxies? If we 
observe neither the GW signal via precision pulsar timing nor inactive close pairs via high angular resolution 
studies, does this imply that some part of the paradigm — the ubiquity of black holes at all redshifts and 
the importance of mergers in galaxy evolution is flawed? Some of these questions have been taken up by 
other groups over the years. In their important paper, Rajagopal & Romani (1995, hereafter RR), discussed 
several mechanisms driving the MBH pair into the GW regime, and calculated the GW spectrum under 
various assumptions. More recently, Gould & Rix (2000) follow up on the gas-dynamical mechanism for 
driving the MBH coalescence first mentioned by Begclman et al. (1980). Armitage & Natarajan (2002) 
explore the detailed interaction of the MBH binary with the accretion disk. Meanwhile, Milosavljevic & 
Merritt (2001) and Yu (2001) have investigated the stellar-dynamical schemes. Finally, Menou et al. (2001) 
have investigated the effect that a time-dependent change in MBH "demographics" might have on the merger 
history of galaxies, and Phinney (2001) has shown an alternative method for calculating the GW spectrum 
for generic sources. 

The advances in our understanding of MBHs along with steady progress in precision pulsar timing have 
led us to this paper. In §2 we first assemble the ingredients needed to model the stochastic GW background 
from MBH binaries, and in §3 we present the results of our calculations. In §4 we discuss the detection of 
gravitational radiation with a Pulsar Timing Array and the status of current experiments. In the concluding 
section we discuss the theoretical and observational prospects for a better understanding of the various 
ingredients to our calculations and measurements. 



First, we need to define the terminology for a stochastic gravitational wave background spectrum that 
has been presented in the literature in a variety of ways {e.g., Burke (1975); Allen & Romano (1999); 
Maggiore (2000)). A single GW is denoted by the transverse-traceless part of the metric perturbation, 
hab{^,t) = /i+(t)e^j,(k) -I- /ix (i)e^j(k), where plus signs and crosses refer to the two polarization degrees of 
freedom of a gravitational wave, and the e^;^ arc basis tensors for the polarizations which are functions of the 
direction of propagation of the gravitational wave, k. The two amplitudes combine with the two coordinate 
directions and the polarization position angle to yield five independent parameters of the single-wave metric 
perturbation. The stochastic background will excite both wave components equally and randomly (|/i+|^) = 
{\hx P) and (/i^/ix) = 0, and all directions and position angles will be equally likely. 

The total power spectral density of gravitational waves, Sh{f), is defined by 



^We refer to this as background radiation, although those interested in relic radiation from the first inflationary moments 
might choose to call this foreground radiation. 



2. The Binary MBH Gravity Wave Spectrum 
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where hp{f) is the Fourier transform of the P = +, x polarization component of the metric strain tensor at 
frequency /, and S is the Dirac delta function (which comes from our Fourier conventions and the definition 
of Sfiif) as a spectral density with units of inverse frequency). There are two spectrum conventions in 
the literature: "two-sided" (defined for |/| < oo) and "one-sided" (/ > 0). They diff'er in amplitude by a 

factor of two: S'onc = 2S'two- Wc use the onc-sidcd version in this paper. From this, wc can define other 
useful quantities, such as the fractional contribution to the energy density of the Universe from a logarithmic 
interval of frequency, 

1 dpGw 27r^ 
7c d\nf " 3HI 

where Hq = 100 ho km s~^ Mpc~^ is Hubble's constant, and pc = SHq/SttG is the critical density in an 
FRW (Freedman-Robertson- Walker) universe. Finally, we can also write down an expression relating the 
power spectrum to a "characteristic strain" spectrum, hc{f), 



ncMf) = -^ = ^fs,m (4) 



hcif) = ^/Tsi^U) . (5) 

We return in §4 to a discussion of the spatial spectrum of the stochastic background as we formulate the 
response of the Pulsar Timing Array to these waves. 

In the following section, we will calculate the number density of sources in a given interval of strain and 
frequency, N{h, /) dh df.^ This is related to the strain power spectrum by 



Sh{f)= / hi^,N{h,f)dh. (6) 
Here, /irms is defined so that the integrated power from a single event counted by N{h, /) is J df Sh{f) ~ 



rms ■ 



Thus, we need to combine various ingredients to calculate the observed gravitational wave spectrum 
due to binary MBHs: the galaxy merger rate, the black hole population demographics amongst galaxies, 
MBH binary dynamics and MBH binary gravitational wave emission, all of which enter N{h,f). In order 
to calculate N{h, /), we first consider a related quantity, N{z) dz, the rate of coalescence events happening 
in redshift interval dz. Since this is a rate, the number of events of a given type in a time bin dt is 
N{z) dz dt = N{h, /) dh df. We can relate this number to what we want by 

.w, ,v/ .dzdt •/ sdz dt dt„ dfr, 

where the subscript p denotes proper values in the rest frame at r{z). Hence, 

N{h,f) = N[zm- (^/.-^ j ^ . N[zih)]-rc^'-^, (8) 
where we have introduced the gravitational wave timescale measured in the rest-frame: 

TGW = (U^) . (9) 



In this derivation we have assumed that the formation rate of the MBH binaries changes slowly on cosmologi- 

cal timescales, and that the binary enters the gravitational wave regime rapidly, so the rate is only a function 
of redshift and not, say, of the initial time of the halo merger. This should be an excellent approximation, 
as argued in Phinney (2001). 



''Here and throughout, wc use the notation Ar(') to refer to the general concept of a number density function, and not to 
some specific function of the arguments. Thus, N(^z, f) dz df refers to the number density of binaries in a redshift and frequency 
interval, whereas N{h, f) dh df refers to the number density of binaries in a strain and frequency interval. 
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2.1. Galaxy Merger /Black Hole Coalescence Rate 

We would therefore like to calculate N{z, Mi, M2) dz dMi dM2, the rate of black hole binary coalescence 
events observed at 2: = occurring in a given redshift interval dz at z, due to MBHs in mass intervals dMi 
and dM2- Because the formation mechanism for these objects in galactic halos is at best poorly understood, 
we will make a simplifying ansatz; 

iV(..M..M,)=.M ^<"-WM'-'' . (10) 

where i'{z) is the total merger rate of the spheroids in which the MBHs reside in a notation similar to that 
of RR and (f^iM, z) is the mass function of black holes at redshift z (normalized so / dM(j){M, z) = (f),{z), 
the number density of MBHs at z). That is. wc assume that the coalescence rate of black hole binaries is 
both independent of their masses and the properties of the halos in which they are present, and rapid on 
cosmological timescales. We will see that we are primarily interested in integrals over the full distribution, in 
which case the individual factors here can be seen as appropriately weighted averages over the multivariate 
distribution. We return to the determination of the black hole mass function in §2.2, and proceed next with 
formulation of the galaxy merger rate. 

We start by considering a shell of some proper depth dip = cdtp at redshift z. This shell has a proper 
volume dVp, or comoving volume dVc of 

dVp = A-Kd\ dip = (1 + z)-^ dVc, (11) 

where dA is the angular diameter distance. This distance is related to the FRW coordinate distance, r, by 
dA = H^^[{aoHor)/{l + z)]; note that, for example, Peebles (1993) refers to r{z) itself as the "angular size 
distance" . The FRW distance is a function of z and cosmological parameters as given by the dimensionless 
integral 

" c r dz' 



aoHoTZ ^ 
aoHor{z) = Sk 



(12) 



_aoHonJo E{z')_ 

where Hq = 100 ho km s~^Mpc~^ is the Hubble constant, ao is the FRW scale factor today, (a/a)^ = 
HqE'^{z) = Hq [flmi'i- + z)^ + Oa + (1 - - f^A)(l + z)"^] , d = da/dtp, Sk{x) is (sin a;, x, sinhx) for (closed, 
flat, open) geometries, and TZ gives the curvature of space. In a flat Universe to which will specialize hereafter, 

aoHoriz) = ^ (flat). (13) 
Jo ^{z ) 

We deflne the merger event rate per comoving unit volume per proper unit time at coordinate time t as 
R{t). The total event rate for the shell is R{t) dVc- The observed rate from the shell requires correction of 
R for the redshift of the time interval, dtp/dt = 1/(1 + z), and then substitution of earlier results: 

g{t) dt = R{t) dVc dtp = R{t) dtp 47rc(l + zfd\ = dtp 4ttc^ Hq"^ {aoHorf . (14) 

Converting the differential from time to redshift, we flnally have 

u{z) dz = g{t) dt = ^^^R{z)H^^\^^^^dz. (15) 

{l + z)E{z) 

That is, J g{t) dt = J v{z) dz gives the rate observed at ^ = 0, with contributions from all along the past 
light cone. This is just the same as equation (12) of RR, except that we write our event rate as R per 
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unit proper time rather than their Fm per unit redshift, although we calculate it as a function of z (so 
Fm dz = R dtp). Also, they specialize to a flat universe with no cosmological constant. In particular, their 
equation (14) is equivalent to our expression with R{z) equal to a constant. 

Carlberg et al. (2000), following the method of Patton et al. (2000), estimate the galaxy merger rate 
out to 2; ~ 1 in the CN0C2 and CFGRS redshift surveys. They count the number of "kinematically close 
pairs" , those likely to be gravitationally infalling, based on their positions and redshifts. They find that 
roughly half of "close" pairs observed in projection at 20h^^ kpc are physically close based on morphological 
disturbances of the pairs. The number of pairs detected is close to what is expected from the two-point 
correlation function. Because the total number of pairs in their sample at this separation is small, they 
determine the overall fraction at this separation by extrapolating the two-point correlation function of 
galaxies from measurements out to 100/i~^ kpc. Using dynamical arguments and results of simulations, 
they find that the average time for a merger starting at this separation of 20h~^ kpc is T = 0.3 Gyr, and 
that a fraction ~ 0.3 of such pairs will merge in this time. Approximately, then, J^/T ~ 1 Gyr~^ gives 
the fractional merger rate for these "close pairs," which should hold roughly constant as long as the galax;y 
population is similar to that of the low-rcdshift universe. 

Using these arguments along with observations of the number of close pairs from the CN0C2 and 
CFGRS surveys over < z < 1, Carlberg et al. (2000) find 



R{z 



b,{z)Rg{z) ~ (j),{z){Qm ± 0.05)(1 + z 



,0.1±0.5 



03 



0.3 Gyr 

T 



Gyr- 



(16) 



where Rg is the merger rate for a single object and 4',{z) is the comoving volume density of merging objects 
which for us is spheroids containing MBHs. We will from here on assume the merger-rate parameters to be 

equal to these canonical values given in equation (16). Crucially, we also assume that the measured merger 
rate per galaxy is equal to the desired merger rate. Inserting this merger rate into equation (15), we find 



v{z) dz = 0.03 yr" 



10-3/ig Mpc" 



0.09 Gyr" 



Ro 



|2 1 



{l + z)E{z) 



dz . 



(17) 



We have taken 0, = 1.0 x 10~^hf, Mpc^'' as our value for the number density of merging objects; see §2.2. 
We have left the evolution of the merger rate, Rg{z)/Ro [with normalization Rp — Rg{0)], unspecified. 
Whereas Carlberg et al. (2000) find Rg{z) approximately constant, Patton et al. (2002), a very similar group 
of authors using somewhat different data, find Rg{z) oc (1 + z)^-^"^-^. To parameterize this difference (as 
well as other substantial observational disagreement on the merger rate as a function of time), we consider 
several other models, generically parameterizing the merger rate per unit time, Rg{z), as a power-law in the 
expansion factor, Rg{z) = Ro{l + z)'*. We will also allow the merger rate per unit redshift to be a power-law, 
which means that Rg{z) = RoE{z){l + z){l + z^'-^l"^, where wc use 7' — 5/2 since that has 7' = 7 when 
= 1 a-iid Q\ = 0. In Figure 1 we show the merger rate with 7 = 0; increasing this exponent would 
increase the merger rate at higher redshift. Below, we will find that our limited knowledge of the merger 
rate is a main contributor to the uncertainty in the final GW spectrum. RR consider merger rates with both 
higher present-day normalizations and stronger redshift evolution (see Figure 1). 
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Fig. 1. — Various scenarios for black hole binary coalescence rates. Light solid line: equation (17) (with 
fim = 0.3, Oa = 0.7, 7 = 2.0); heavy solid line: equation (17) {Qm — 1-0, JIa = 0.0, 7 = 2.0); dash-dotted: 
equation (12) from RR; dashed: equation (13) from RR. 
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We note as an aside an alternative to those' phcnomcnological and observational merger rates. Menou 
et al. (2001) have instead considered theoretical calculations of the merger rate of galaxy halos. They perform 
Monte Carlo simulations of the merger history of halos using semi-analytic methods (so-called "merger trees") 
based on numerical simulations and the Press-Schechter formalism and its extensions. This allows them to 
vary the black hole mass function and its relationship to the underlying halo mass function as a function of 
time. In this work, however, they do not consider the detailed gravitational wave spectrum that results from 
their models. 

Our model is useful beyond the nHz-/xHz regime. The LISA satellite will be sensitive to the total event 
rate of MBH binary inspirals as they go through their final coalescence. This event rate is just J ^{z) dz, 
with i'{z) given by equation (15) or equation (17). Wc show this quantity in Figure 2. For the models we 
have been considering, the event rate is 0.01-1 binary coalescence events per year. However, only a subset 
of these events will be detected with sufficiently high signal-to-noise for a sufficient length of time over the 
different phases of the MBH binary merger to enable a measurement of their individual masses and distances 
(e.g., Hughes 2002). 
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Fig. 2. — Total MBH binary coalescence event rate, for two sets of cosmological parameters and two choices 
of parameterizing galaxy merger rate, parameterized by different maximum redshifts and power-law indices 
7 (black, solid) and 7' (red, dashed) as marked. (In the right-hand panels, the curves lie atop one another.). 
We have assumed a present-day merger rate of i?o = 0.09 Gyr^^ and MBH density of = 1G~^/iqMpc~^. 
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2.2. Black Hole Population Demographics 

Magorrian et al. (1998) showed that essentially all "spheroids" (elliptical galaxies or the bulges of spirals) 
of mass Afgph have central black of mass M, ^ 0.006A'/sph- More recently Gebhardt et al. (2000), Ferrarese 
& Mcrritt (2000), Merritt & Ferrarese (2001a) and Tremaine et al. (2002) have shown that the data are 
consistent with an even tighter correspondence between the black hole mass and the velocity dispersion 
of the galaxy raised to a power between about 4 and 5. These authors find a significantly smaller ratio of 
M, to Mgph than that in Magorrian et al. (1998) with more careful treatment of the data. The a^, — M, 
relation further implies that MBHs may grow primarily by other process such as accretion, rather than by 
the coalescence process discussed here (Merritt & Ferrarese 2001b): if the spheroids grow simply by mergers, 
and the MBH mass initially tracks the spheroid mass, we might expect that the MBH mass would grow 
linearly with the mass of the spheroid, whereas the velocity dispersion raised to some power would not. 
There is some controversy, which is discussed by Merritt & Ferrarese (2001a), of the implications of these 
results for the detailed M.-Mgph relationship. 

Determination of either the mass function, the velocity dispersion function, or the luminosity function 
and the mass-to-light ratio for the population of galaxies with spheroids will provide us with the mass 
function of central black holes. Unfortunately, there have been no good measurements of the velocity 
dispersion distribution for galaxy spheroids. Instead, we use the original, less tight, correlation between 
black hole mass and spheroid mass. We will take the differential probability of black hole mass at a given 
spheroid mass, (j){M,) dM,, to be lognormal with mean and dispersion given by 

(logio(M./M,ph)) = logio(0.0012) ± 0.45 (18) 

(Merritt & Ferrarese 2001b). We then combine this with a theoretical model for the spheroid luminosity 
function (Fukugita et al. 1998; Magorrian & Tremaine 1999, Maggorian, private communication) which is 
also lognormal with 

(logio(Lsph/i0)) = logio(2 X 10*^) ± 0.6 . (19) 

Because we want an actual spatial density, we must normalize the luminosity function with the spheroid den- 
sity, (j), = 10-3 hi Mpc-3. We also need the spheroid mass-to-light ratio, M/L = 7.0 {L/10^°Lq)° '^Mq/Lq. 
Putting all this together by integrating over the two lognormal distributions gives us the MBH mass function, 
also in lognormal form, with 

(logio(M.Af0)) = logio(1.2 X 10^) ± 0.6 (20) 

This expression is appropriate for elliptical and SO galaxies, and for the spheroidal component of disks and 
spirals. There is also evidence {e.g., McLure & Dunlop 2001, 2002) that the same or very similar relations 
hold for low-redshift {z < 0.2) AGNs (Seyferts and QSOs), although there is some dispute regarding the 
constants of proportionality. 

In Figure 3 we compare this mass function with the earlier model of Small & Blandford (1992), which 
was used in RR. Although the shapes of the distributions are similar, the more recent results on demographics 
predict that MBHs are more common: there are roughly an order of magnitude more at masses below IO^Mq. 

This hard data on central black hole populations is of necessity at low redshift. In this work we are also 
concerned with the evolution of this population. The very existence of active galaxy populations at high 
redshift of course implies that some fraction of galaxies have had central black holes for a long time, but 
the detailed demographics are less certain. In their recent work, Menou et al. (2001) consider a scenario in 
which the fraction of spheroids containing an MBH is a function of redshift. In the following, such evolution 
is completely degenerate with changes in the merger rate, but for definiteness we will consider two scenarios. 
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one in which the black hole mass function is constant with time, and another in which the normalizing 
number density follows a power law in {1 + z); this will allow us to combine a possible evolution in the 
merger rate with evolution in MBH populations as a single power law as in the previous subsection. In 
general we would certainly expect the peak of the mass function to also evolve over cosmological timescales, 
an issue we defer to later study. 
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Mass [M/Mq] 



Fig. 3. — The black hole mass function; the smooth upper curve is the lognormal model presented in the 
text; the black jagged lower curve is from Small & Blandford (1992), as used in Rajagopal & Romani (1995). 
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2.3. Black Hole Binary Dynamics 

Once the merger of two galaxy halos occurs, we must then follow the evolution of the central black holes. 
This is the subject of some debate. The paucity of obvious close binary nuclei imply that mergers between 
black-hole possessing-galaxies must be rare or short-lived. Wc know that the host galaxies themselves have 
often experienced at least one "major merger" in their lifetime. This is especially true of the most massive 
ellipticals harboring the most massive black holes. The recent study of the dynamics of MBHs in merging 
galaxies given by Yu (2001) builds on the seminal work of Begelman et al. (1980). 

The evolution of the MBH binary proceeds in three stages. First, the pair spiral into what becomes the 
common single core of the merger remnant driven by dynamical friction, on a timescale given by 

„ / 4xl06yr \ / U 

logW* y V200 kms-ini00pcy \10H'IqJ ' ^ ' 

where the core of the galaxy has radius Vc, velocity dispersion Uc and contains iV* stars, and M, is the mass 
of the smaller black hole. 

The binary continues to harden (inspiraling and gaining energy) first through dynamical friction, sub- 
sequently by three-body interactions, and finally as a hard binary, when the semimajor axis is given by 

and the Keplerian period is given by 

^^^•^>^10'^^(l0^)(200^)"- (23) 

At this separation, however, the timescale for evolution via gravitational radiation, given by equation (2), 
is much longer than a Hubble time. So the typical pair must somehow evolve to a much tighter orbit with 
a ~ 0.02 pc in order to ever complete the merger. At this separation the Keplerian period is P ~ 30 yr 
and orbital speed is -^4000 km s~^. Various dynamical mechanisms have been proposed for this process. 
Begelman et al. (1980) first discussed the dynamics in detail and noted the difficulty of reaching the rapid 
GW inspiral regime, as the MBH binary empties the "loss cone" in the stellar distribution function of 
stars with which the MBH binary can interact. RR surveyed the field at the time, and decided that while 
"unassisted stellar dynamics" would bring only a small fraction to the GW regime, "external perturbations 
may, however, cause efficient inspiral." Examining these possible perturbations further, Gould & Rix (2000) 
and Armitage & Natarajan (2002) revisited the earlier idea (Begelman et al. 1980) that gas dynamics — the 
same physics as planetary migrations — could be responsible for this evolution, and indeed for fueling some 
or all quasar activity. Recently, Milosavljevic & Merritt (2001) have done detailed N-body calculations of a 
MBH-MBH binary falling into a spherical potential of stars, and discuss the issue of potential stalling of the 
binary coalescence. Yu (2001) has examined in yet greater detail the stellar-dynamical schemes in tri-axial 
galaxies. Yu concludes that indeed there should be a sizable population of binary MBHs remaining in massive 
post-merger galaxies, in the abscnicc; of such gas-dynamical effects, but notes the difficulties of detecting the 
remnant binaries. On the other hand, Zhao et al. (2002) argue that the central stellar distribution of galaxies 
is considerably more complicated (cuspy and asymmetric) keeping the loss cone populated, possibly rescuing 
the stellar-dynamic scenarios. 

Another possibility is that the binary will be driven to the inspiral regime by the presence of a third 
MBH from a second merger event. Tidal forces from the third body — or equivalent perturbations from a non- 
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axisymmetric galactic potential — can decrease the merger time by an order of magnitude and significantly 
increase the eccentricity of the binary via the "Kozai mechanism" (Blaes et al. 2002). 

In the following, wc will assume that some mechanism like gas physics is indeed successful. Thus, all 
relevant pairs enter the gravitational wave regime and eventually shed their orbital energy in the form of 
gravitational radiation. The final ten thousand seconds or so of this evolution are thought to be the most 
energetic GW events in the (present-day) Universe, and will be easily detected by the LISA satellite over a 
wide range of MBH masses and redshifts (Hughes 2002). Here, though, we are concerned with the quiescent 
quasi-Keplerian evolution that precedes the infall events. 



2.4. Gravity Waves from Binary MBHs 

Now, we will briefly review the emission of gravitational radiation from binary black holes. First, we 
will need the amplitude of the gravitational radiation emitted by a binary at distance D. This is given by 
(Peters & Matthews 1963; Thorne 1989) 

where M = [MiM2(Mi + Ma)-^/^]^/^ is the "chirp mass" of the system, and Pp is the proper rest-frame 
period of the binary. The observed frequency of the radiation in the nth harmonic will be / = n/p/(l -|- z), 
where fp = 1/Pp- Only terms with n > 2 are nonzero, and n = 2 for the circular orbits to which we will 
restrict ourselves hereafter (Quinlan 1996; Quinlan & Hernquist 1997). However, if the binaries are driven 
together by the presence of a third MBH (Blaes et al. 2002), the eccentricity may be significant and populate 
higher-frequency harmonics. This would only serve to increase the amplitude of gravitational radiation and 
populate the higher harmonics (by an eccentricity-dependent factor). This amplification would offset, or 
exceed, the loss resulting from expulsion of the third body which is assumed to coalescence in our model. 
Equation 24 differs by a factor of •\/3/4 from that given in RR. Instead of the "characteristic strain" from 
Thorne (1989), which includes this factor to recover the signal-to-noise for an earth-bound detector, we just 
use the angle-averaged mean-square strain [i.e., eq. [24] is l/\/2 times the maximum rms strain). We point 
out for ease of interpretation of this and later formulae that {GM/c^) has units of time. In a cosmological 
setting, we replace D with caoHor{z) / Hq. 

We will also need the characteristic timescale of the emission, defined in equation (9). Using the 
Kepler formula, /^a^ = {2-k)~'^G{M\ + M2) and the formula for energy loss from gravitational waves, the 
gravitational wave timescale is {e.g., Peters & Matthews 1963; Shapiro & Teukolsky 1983) 



_ , 1 96 (GM-\'/\ 



z)f'\ (25) 



where we have used fp = f{l + z)/n and set n = 2 for circular orbits. 

3. Strain Spectrum 

Now we can combine all of the ingredients into an observable strain power spectrum. We have 
• i'{z) dz, the volume merger rate observed at ^; = of galaxies with MBHs located between z and z + dz; 
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• (l>{M,) dM, the distribution of MBH masses, normalized to the number density of galaxies whose 
merger rate is given by iy'{z); 

• hrt^s{z, f, Ml, M2), the strain observed at 2; = for a binary of given masses with frequency / at 
redshift z; and 

• TGW, the gravitational radiation timescale. 

Putting all this together using equation (7), the strain distribution is 

dz 

N{h, /, Ml, M2) dh df dMi dM2 = N{z, f, Mi, M2) —dh df dMi dM2 (26) 
where the distribution in redshift is 



-5/3 



5 Anc^ <t>{Mi)<t>{M2) fGMY'^\ [aoHor{z)f g/a 1 

= 96-Hr ^1 [—) ^^'^E{z)(l + zf/^ ^""^^ /■ (2^) 

First, note that this expression differs from equation (20) of RR; they neglected to include the factor 
of (1 + z) needed to convert from proper rest-frame time interval tp to observed interval t. Second, note 
that equation (27) factors into functions of mass, redshift and frequency. This means that marginal densities 
(integrals over one or more of the variables) are easy to calculate and have the same dependence on the 
remaining parameters as the unmarginalized distribution. This will make the following calculations particu- 
larly simple. However, it is certain that a fuller treatment of the model, incorporating the likely dependence 
of the halo merger and MBH coalescence rate simultaneously upon mass and redshift, including a complete 
accounting for the probability of a given MBH binary of successfully inspiraling to reach the GW regime, 
will no longer display this simplicity. 

We now formulate the strain spectrum using the definitions of equations (5-6): 

Kif? =f Jdh dMi dM2 hl^,{z, /, Ml, M2)N{h, /, Ml, M2) . (28) 

We can simplify this expression considerably by using N{h, . . .) dh = N{z, . . .) dz, which allows us to remove 
the dz/dh factor in equation (26). Thus 

hciff = fjdz dMi dM2 h^^,{z, f, Mi,M2)N{z, f. Mi, M2). (29) 

Using the results of equations (24-29) and the developments of the previous sections, the strain spectrum is 
given by 

K{ff = jdJ-^ [nfil + z)]-'^' iiGM/c^f\j^y (30) 
where we use (• • •), to denote averages over the MBH population. In equation (30) 

(M^f^), =1 JdMi dM /^^^y^^^ A^^/^ = (2.3 X lO'M^f , (31) 
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where wc have used the MBH mass function of §2.2, and have assumed that the number density of MBHs 
and the hales whose merger rate is described by equation (17) arc the same. 

We can now also separate out the z integration to see the simple form of the result: 



R{z) dz 

HoE{z){l + z)*/3 

5/3 / -f \ -4/3 



In the last line we have set the number density of merging galaxies and their present day merger rate to 

(j), = 1.0 X IQ-^hl Mpc"^ , and 

Rg = 0.09 Gyr"^ . (33) 

Phinney (2001) derived an expression equivalent to equation (32), and noted the weak dependence of hdf) 
on cosmology owing to a cancellation of aoHor{z) factors. Only the E{z) factor remains in the integral. 
However, the details of the distribution N{h, f) do depend more strongly on the cosmology. 

In Figure 4 we show the dimensionless integral, 

_ r R{z) dz 

^''=J -R^E{z){l + zr/^' ^^^^ 

which is essentially the same as ((1 + z)^^/'"*) in Phinney (2001). Wc plot 1^2 for a variety of cosmologies, 
merger rates and black-hole mass-function histories. We parameterize the latter two by the single exponent 
7, defined such that (j){Mi, z)(j){M2, z)R{z) 0(.{1 + z^ or oc E{z){l + z)'^'-^/^ as in §§2.1-2.2. The maximum 
redshifts of 3 and 5 used in Figure 4 are motivated by the distribution of objects with redshift which is 
discussed below in §3.2. For 7 < 2, we are most sensitive to z < 3 and the cutoff is irrelevant. For stronger 
merger-rate evolution, we are sensitive to more distant events. Figure 4 shows the weak dependence upon 
both the cosmology {flm and IIa) and the maximum redshift of MBH binary formation. The results also 
show the weak dependence on whether we consider the merger rate as a power law in redshift (with exponent 
7) or time (7'). Similarly, if wc allow this exponent to vary by as much as one or two units, wc again get 
very little change. The dependence on the maximum redshift is obviously somewhat stronger if we increase 
the power-law exponent 7, thereby putting more events at higher redshift. There remains a prefactor of Rq 
in h'^{f) — the overall normalization of the merger rate — which can have a strong effect. Of course, each of 
these effects can change things by a factor of two. In summary, the uncertainty in h"^ (/) is at least an order 
of magnitude. 
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Fig. 4. — Redshift integral occurring in the MBH gravitational wave strain spectrum, equation (34), as in 
Figure 2. The strain spectrum is proportional to this quantity. 
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3.1. The Strain Distribution 

We would like to estimate the entire probability distribution of h'^if) about this mean. The simplest way 
to approach this is to use the method of characteristic functions, which are simply the Fourier transforms 
of probability densities. In order to do this, first consider the quantity 

hl{f)j = J dz dMi dM2 hl^,{z, f. Ml, M2)N{z, /, Mi,M2)df . (35) 

In this equation, the quantity N{z, f. Mi, M2) dz df dMi dM2 gives the actual number in a cell of size 
dz X df X dMi X dM2 at a given frequency and bandwidth / and df. It can be shown (see Appendix A) that 
the characteristic function of the distribution of h'^if) df/f is given by 

liiifit) = J dz dMi dM2 N{z, f, Ml, M2) [e*t''™(^./.Mi,M2) _ -^j ^gg^ 

While one can show that this integral exists for the cases we are considering, it cannot be expressed in 
closed form. For this reason, the distribution itself — the inverse Fourier transform of 4>(t) — is quite difficult 
to calculate. 

To start, we can use the fact that the coefficients of the Taylor expansion oi\n(j){t) around t = Q give Kfe, 
the 'semi-invariant moments' or 'cumulants' of the distribution {i.e., the mean, variance, skewness, kurtosis, 
etc. for fc = 1, 2, 3, 4, . . .), which we can read off as 

Atfe(/) = j dz dMi dM2 [h'i^,{z, /, Ml, M2)] N{z, f. Ml, M2) df . (37) 

Note that these arc the moments of h^{f) df/f, not h^{f) itself, which can only be measured with a finite 
bandwidth. For example, the variance of h1{f )df is 



var 



[hlif) df] = f^ var [hUf) df/f] = f^ J dz dM^ dM2 hi^,N{z, f. Mi, M2) df . (38) 



We measure h1{f) itself with some finite bandwidth. A/, as h1{f) ~ /^^ h^Af) df / ^f . The variance on this 
quantity decreases as the bandwidth increases. 

The variance of the MBH density N{h, f) under consideration is 

var [hl{f)] = Anc^Ho^ (^^) ' {{GM/c'f). J dz [aoHor{z)]-^ R{z)/E{z) . (39) 

As before, this separates into an average over the MBH mass function and an integral over redshift. Consid- 
ering the latter, we see that the integral in question diverges as 2; — > 0, with the integrand behaving as dz/z^. 
Indeed, all of the integrals with fc > 1 diverge: the moments do not exist. Nonetheless, the characteristic 
function and the distribution itself remain well-defined. This divergence occurs because the strain from a 
single event is proportional to the inverse of the distance; although the mean-square accumulated strain (the 
power spectrum) converges, the dispersion from small distances diverges. 

We have performed numerical experiments with the simpler problem of a uniform MBH binary distri- 
bution in a non-expanding, Euclidean universe, so N{z) oc and /!^jjjg(z) cx 1/z^ for all z < z^i^^ — an 
approximation to the situation at z <C 1. We find that the distribution of strains is highly skewed, with a 
power-law tail toward larger amplitudes. This tail is due to the possible contribution of low-redshift pairs 
which give a contribution h'^iz) oc l/z"^. 
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3.2. Other Quantities 

Our derivation also allows us to calculate other quantities related to the distribution of strains. First, 
we can ask what is the number of binaries in a frequency interval, 

N{f)df = J dzdMidM2N{zJ,Mi,M2)df . (40) 

Using equation (27), and again under the assumptions we have made so far that the (f), is independent of 
redshift, the z integration separates out, and we see first that N{f) df/f oc f'^'^df/f, a steeply-falling 
function of frequency. In full, 

M y"\ ( f Y'^' [ ^ R{z) [a^H^r{z)f 



Ar(/)./ = 5.6xlO-W(^ ^^^^^,^J j'^^ Eiz^ilizr^ - ^''^ 

{M.-^'^), is (4.1 X l{fMQ)~^'^ for our MBH mass function which is different from the value above owing 
to a different weighting function. In Figure 5 we show the dimensionless integral occurring in this expression 
for a variety of cosmologies and merger histories. The stochastic GW background at nHz (yr~^) frequencies 
is therefore the result of nearly a million simultaneous binaries across the Universe. 
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Fig. 5. — The dimensionless integral occurring in N(f), 
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jquation (41), for a variety of scenarios as in Figure 2. 



- 21 - 



Similarly, we can ask, what is the range in redshift that we are probing? This is just 

i?o E{z){l + zf/^- ^ ' 

which wo show in Figure 6. For larger values of the exponent (7 > 2), relatively more events occur at 
higher redshift, and the dependence of the observed strain on the maximum redshift will be stronger. More 
specifically, we can calculate the mean redshift, 

^'^^^^ - SdzNizJ) ■ (^^^ 

Again, the integrals separate and we find that this average redshift is independent of frequency: 

^ Jdz R{z) [jaoHorf/Eiz)] (1 + z)-^/^ _ 

JdzR{z)[{aoHor)yEiz)]{l + z)-»/3 ' ^ ' 

This is independent of the details of the MBH population and the normalization of the merger rate, but 
dependent on the merger rate evolution and the cosmology In Figure 7 we show the average redshift for a 
variety of cosmologies and merger histories. 
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n^=0.30, 0^=0.70, 7 = 0.0 




z 



n^ = 0.30, 0^=0.70, 7 = 2.0 




z 



Fig. 6. — The distribution of redshifts probed by MBH Binaries, equation (42), for two sets of cosmological 
parameters and two choices of parameterizing galaxy merger rate. The sohd Hues parameterize the merger 
rate as a power-law in redshift; the dashed lines as a power law in time, as in Figure 2. 
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Fig. 7. — The average redshift from which the MBH gravitational wave signal comes, equation (44), for a 
variety of scenarios as in Figure 2. 
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In the previous sections we calculated the strain spectrum, h'^{f), and variants. Using the above 
developments, we can calculate a slightly different quantity: the average strain amplitude that contributes 
to hlif). This is just 

{hV)) = §^^f'^'' (45) 

which rises owing to the more rapid decrease of number density with frequency relative to the spectrum 
amplitude decrease. That is, higher frequency signals are from small numbers of rare, powerful events in 
comparison with lower frequency signals. 

Finally, we can also calculate quantities related to the mass distribution that contribute to the strain 
spectrum. Again due to the factorization of N(M\, M2, z, f) these are independent of redshift and frequency. 
For any function of the masses, we see that 

where as above (•••). refers to averages over the MBH distribution. One particular quantity of interest is 
simply the average chirp mass: 

{M) = 2.9 X 10^ Mq. (47) 

which, as would be expected, is close to the average of the MBH distribution itself (Fig. 3). Another is the 
mass ratio, 

(q) = (M./M^) = 11 . (48) 

Because the MBH mass distribution is itself so wide (Fig. 3), the gravitational waves sample a fairly wide 
range of black-hole mass ratios. In particular, we are sensitive to minor mergers as well as major mergers. 
However, these quantities change to {A4) = 5 x lO^M© and (q) = 1.7 if a minimum black hole mass of lO^M© 
is assumed. Such a minimum hole mass is warranted based on considerations of the strong dependence of 
dynamical friction time scale on mass (Xu & Ostriker 1994). The values in equations (47-48) can be expressed 
as individual masses (Mi, M2) ~ (1.3 x lO'^, 1.1 x 10*^) Mq, or (7.4 x 10^ 4.4 x 10^)Mq if the distribution is 
cut off at 10^ Mq. 

In a similar vein, we can weight these functions by the signal amplitude they produce, rather than the 
number of binaries contributing. Since /i^^g oc M.^^^^ (eq. [24]), this changes the exponents in equation 46 
from —5/3 to +5/3. Such a weighting gives a signal- weighted average chirp mass of 9 x 10 ""M© and mass 
ratio of 6. Again, we see that a small number of high-mass black holes are responsible for the bulk of the 
signal. Similarly, the signal-weighted average redshift is somewhat lower than the number-weighted redshift 
of Figure 7. 



3.3. Results 
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Fig. 8. — Characteristic strain spectrum hc{f) for the fiducial models discussed in the text, along with 
Monte-Carlo realizations. The upper thick curve and the associated realizations have 7 = 2; the lower set of 
curves have 7 = 0. The blue dashed line gives the current best limits on the gravitational wave background 
from pulsar timing observations. The dotted line shows the expected limits from a Pulsar Timing Array, 
after operation for ~ 8 years. 
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In Figure 8 we consider two models of galaxy mergers and the black hole population, and show the 
resultant spectrum of characteristic strain along with several realizations of a Monte-Carlo simulation of the 
spectrum. We perform the Monte Carlo by simply Poisson-sampling from N{z, /, Mi, M2) on a grid of the 
parameters and summing the contributions of h'^^g{z, f, Mi, M2). We have chosen flrn = 0.3, Oa = 0.7, the 
constant halo merger rate given in equation (17), the MBH mass function from equation (20), and a maximum 
redshift of Zmax = 3. The upper set of curves have a relatively strong power-law index 7 = 2 describing the 
galaxy merger rate; the lower curves have 7 = 0. We have assumed a bandwidth of Sf/f ~ 0.12. We note 
that at higher frequencies the simulation is somewhat below the expected mean spectrum. This is due to 
the high-amplitude tails mentioned in the discussion of the full strain distribution, §3.1 above: small number 
statistics and numerical inaccuracies mean that we miss the large spikes. 

We expect the actual dispersion to be greater at higher frequencies — smaller numbers of objects are 
contributing more due to the steeply-falling function N{f) (eq. [40]). However, just as the mean of our 
simulation falls below the expectation value, so does the dispersion around the mean: we miss the large 
excursions. This is a larger effect for the relatively high value 7 = 2 than at 7 = 0, representing comparatively 

fewer merger events at low redshift. 

However, a central result of this paper is that the slope of the spectrum is independent of all of this. 
As in Phinney (2001), we find that hUf) cx or Shif) oc /^^Z^, or nawif) oc p/^. The overall 

amplitude, however, will vary up and down by the factors calculated above and shown in Figure 4. As 
mentioned above, this means that, on the one hand, we cannot use the information from any observed slope 
to limit the physics. On the other hand, it potentially provides a check that we are observing the background 
from a coalescing population (although it is not sensitive to the details of that population). Indeed, other 
proposed sources of stochastic gravitational waves often have different slopes. Early universe cosmological 
backgrounds, for example, generically have fi(/) oc const and hence h^{f) oc /~^. 

4. Gravity Wave Detection with a Pulsar Timing Array 

We turn now to the experimental approaches to detection of the stochastic background from MBH bina- 
ries. We begin with a statement of the response of pulsar timing to a plane GW. In the following subsection 

we discuss the angular and temporal basis functions required for detection of a stochastic background. In 
the final subsection the current limits and future prospects are discussed in relation to the expected signal 
level presented in Figure 8. 

4.1. Plane Wave Response 

Precision timing of the arrival of pulses from a pulsar allows the direct detection of gravitational radiation 
traversing the sightline (Sazhin 1978; Detweiler 1979). The spacetime metric strain perturbation h can be 
viewed as a perturbation from unity of the index of refraction for electromagnetic radiation in flat space. 
For a plane GW traversing the entire sightline to a pulsar the pulse arrival time is advanced, or retarded, 
in proportion to the incomplete cycles of the wave traversed by the pulse at the pulsar and in the solar 
system. Observable effects are possible only when the GW amplitude, or state of the spacetime metric, 
changes locally or at the pulsar or both in ways distinguishable from pulsar spin and other parameters. A 
useful analytic approach is to formulate the apparent redshift Z of a pulsar spin frequency which arise from 
gravitational radiation. Observationally Z is estimated by taking the derivative of the residuals that result 
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after fitting pulse arrival times to a model for the pulsar's spin properties and astrometric coordinates, and 
binary parameters if relevant. For a single plane wave the observable may be written as 

Z{a,0,r,tr) = ^'^.~fh h+{ctr) - h+{ct^ - 70] + --^^[/ixK) - h^{ct^ - 7O], (49) 
2(1+7) (1 + 7) 

where (a,/3, 7) are the direction cosines of the pulsar with respect to the GW, {h+, hx) are the polarization 
components of the GW, and {te,tr) are the times of emission (e) and reception (r) of the pulses from the 
pulsar which is at a distance /. 

The ti terms in equation (49) will be correlated between different pulsars. Therefore timing measure- 
ments of an array of pulsars across the sky, the Pulsar Timing Array, creates a gravitational wave telescope 
that is sensitive to a spectrum of waves with periods less than the duration of the measurements which is 
years, or 10 nHz {e.g. Backer 1993). The five degrees of freedom in equation (49) correspond to the five 
degrees of freedom in the trace-free space-time metric. The terms in equation (49) will be uncorrelated 
between different pulsars, and yet will have comparable amplitude. These terms create an irrcdiicible noise 
background in addition to any limitations from the measurement errors and properties of the pulsars them- 
selves or the intervening turbulent plasma through which the signals propagate. Lommen & Backer (2001) 
apply these ideas to a search for plane waves from both our Galactic Center and nearby galaxies which host 
MBHs based on the hypothesis that the central objects are binary. 



4.2. Fitting Algorithms for a Stochastic Background 

In the preceding sections we considered the production of a stochastic background of gravitational 
radiation from coalescences of MBH binaries. This is summarized in the form of a characteristic strain 
spectrum. h,.if), Figure 8. For the purposes of analysis of timing data from a spatial array of pulsars 
we need to consider the full tensor perturbation field hij (x, t) . All elements are statistically independent 
owing to the random superposition of many sources from many directions. This superposition destroys the 
simple angular plane wave pattern described in equation (49). However, Burke (1975) has shown that 5/8 of 
the variance in the stochastic fluctuations can be extracted using the 5 quadrupole spherical harmonic Yj^ 
functions as the angular basis. 

Temporal modulations of each spherical harmonic term can be described using either Fourier frequency 
coefficients as the basis, or polynomials (Foster & Backer 1990), or orthogonal polynomials (Stinebring et al. 
1990). The polynomial approaches are particularly suited to this analysis owing to both the steep spectrum 
expected for the MBH-MBH GW emission (eq. [32]) and the need to fit for the spin parameters of the stars. 
Given any choice of the temporal basis function, the amplitude /irins(/) for any temporal term is formed 
by the square root of the sum of the squares of the 5 spherical harmonic coefficients. Conversion to hc{f) 
requires consideration of the spectral window function which will be addressed in a later work. 

A second approach to the angular basis has been suggested by Hellings (1990). The tensor field stated 
above can be transformed into a spectrum of waves /iij(Q, /). These can be decomposed into waves traveling 
along the three cardinal directions, each with independent polarizations: 

hii hi2 hi3 
hi2 /122 /l23 
. /ll3 /l23 

In this case one of the six terms in the decomposed format is not independent owing to the trace- free property. 
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The Doppler patterns of the three waves can then be summed to create Doppler patterns for each tensor 
element. Thus the full variance of the stochastic background can be sampled. As in the spherical harmonic 
case each term will be temporally modulated, and the quantity /irms is formed by the square root of the sum 
of the squares of 5 independent terms. 

4.3. Current Limits and Future Prospects 

While the preceding subsection outlines approaches to detection of the GW background using an array 
of pulsars, measurements from a single object can establish an upper limit on the background. Kaspi et al. 
(1994) has stated the best limit at nHz frequencies using precision timing of PSR B1855+09 relative to the 
world's best atomic time scale. A second pulsar in their study, PSR B1937+21, was shown to be unstable 
on long time scales, an effect they attribute to internal structure of that neutron star. Thorsett & Dewey 
(1996) improved the limit using further statistical considerations. At a typical frequency of 5 nHz (8 yr 
period) they conclude that the energy density in gravitational radiation per logarithmic frequency interval 
is less than Q,gh? 6 x 10~^, where h in this expression is the ratio of Hubble's constant to 100 km s~^. 

Energy density scales quadratically with the time derivative of the metric strain. Therefore the level of 
ilg/i^ that can be detected by pulsar timing scales quadratically with the rms timing residual i?, and as 
with the duration of the experiment T, which is the inverse of the minimum frequency sampled. Recently 
Lommen & Backer (2002) have extended the Kaspi study of PSR B1855+09 by more than doubling the 
experiment duration to 17 y. The new upper limit on energy density is roughly an order of magnitude lower 
at frequencies near 5 nHz. A full report on this work is in preparation. An approximate statement of the 
Kaspi-Lommen limit on characteristic strain is given in Figure 8, which shows that the measurements are 
approaching a level of significance given our current model of the Universe. An alternate way of stating this 
is that the measurements are placing useful constraints on the uncertain parameters in our model. The ratio 
of the pulsar detection level to the characteristic strain spectrum level scales as Ti3/6_ 

The Pulsar Timing Array experiment, which will use both existing and new data sets, can improve on 
the Kaspi et al. (1994) result in three obvious ways: (a) smaller timing residuals owing to combined data sets, 
new and upgraded telescopes and better data acquisition techniques; (b) more objects which both provide 
the capability of actual detection as opposed to upper limits and, if sufficiently numerous, provide a "root 
N" advantage; and (c) longer experiment duration. In Figure 8 we show a reasonable goal for the future: 200 
ns timing precision over eight years. Current measurement series are already achieving this level of precision 
for a few objects. 

5. Summary 

We have calculated the spectrum of the stochastic background of gravitational radiation from the Uni- 
verse of coalescing binary black holes in the centers of galaxies with simple parameterizations of the current 
uncertainties. This is followed by a discussion of the influence of the stochastic background on precision tim- 
ing measurements of pulsars which includes upper limits on the background and use of an array of pulsars 
for direct detection. This work was motivated by improvements both in our knowledge of the present-day 
massive black hole population and the rate of galaxy mergers as well as new pulsar measurements. 

Following the approach of past authors, we construct the spectrum from coalescing binary MBHs from: 
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the galaxy merger rate, the black hole population demographics amongst galaxies, MBH binary dynamics 
and MBH binary gravitational-wave emission. Our galaxy merger rate is based on observations of close 
pairs and an estimate of their dynamical friction time scale. We convert from the galaxy mass function to 
the black hole mass function using recent determinations of the correlation between black hole mass and 
spheroid mass. 

Our "fiducial" model has rapid evolution in the merger rate per unit time per galaxy as a function of 
redshift, and a constant mass function of MBHs out to z = 3. The characteristic strain spectrum predicted 
for this model, hdf) ~ 10^^®(//yr^^)^^/'^, is just below the latest observational hmits at / ~ 0.2 yr^^. 
Whereas the slope and general character of this prediction agree with the earlier work of Rajagopal & Romani 
(1995), our predicted amplitude is somewhat higher. This increase is the result of the order of magnitude 
increase in the local MBH number density, which is tempered by a somewhat lower merger rate at high 
redshift. 

With our formulation we can calculate other quantities. The number of binaries contributing for unit 
bandwidth {Sf/f^ 1) at nHz frequencies is approximately 10^. It is evident from our simulations that the 
variance in the strain spectrum in the same frequency band is roughly 50%, considerably larger than the 
10"'^ that would be expected if all events were weighted equally. This indicates that the spread around the 
mean is due to the possibility of a small number of high-amplitude, nearby events. The median redshift for 
the binaries contributing to the spectrum is < 1 if the merger rate evolves relatively slowly with redshift 
(7 < 2), but can be large for strong evolution of the merger rate. The mean "chirp" mass is 3 x 10^ M© and 
the mean mass ratio is 11; the individual masses arc then 10^ and 10*' Mq. 

Our calculations allow us to see directly how the remaining lack of understanding of some of these 
physical processes impacts our predictions. These are areas for future research. The most important factors 
are: 

• The value of the present day galaxy merger rate is a matter of considerable debate. The value we have 
quoted, 0.09 Gyr~^ per galaxy from Carlberg et al. (2000) is known to no better than 50%. Moreover, 
this value is measured for massive L* galaxies; does it apply to the wider range harboring MBHs? 

• The evolution of the merger rate for moderate redshifts (^ ^ 1) is even less well-known. Looking back 
still further, the evolution of the merger rate at z > 1 is yet more difficult to pin down. The high 
redshift surveys, in progress and planned, arc crucial for progress. 

• We have not taken into account the details of MBH binaries dynamics: under what circumstances, if 
any, do MBH binaries make it to the GW regime? We have assumed that all galaxy mergers promptly 
lead to coalescence of their central black holes. A more complex model needs to consider both a delay 
between merger and coalescence which may depend on mass and the possibility of expulsion in a triple 
system. 

• Throughout, we have simplified the MBH coalescence rate by splitting the MBH population from the 
galaxy merger rate, as in equation (10). Although this ansatz is useful for examining the dependences 
of the results upon the various ingredients, a full treatment of the time- and mass-dependent merger 
rate of galaxies (or the black holes directly) needs to be included. Appropriate results could come from 
considerably improved observations, N-body simulations, or formalisms such as the Press-Schechter 
calculation of the halo mass function. 

• Finally, we have assumed a simple model for the black hole population, combining the spheroid mass 
function and MBH demographics, as in §2.2. This does not take into account more recent results 
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relating the MBH mass with the gravitational potential of its parent spheroid via the velocity dispersion 
(Gebhardt et al. 2000; Ferrarese & Merritt 2000). Moreover, we assume that the evolution in the MBH 
mass function can be modeled simply as a scaling of the number density, rather than the perhaps more 
realistic scenario where the mean and spread of the mass also evolves. 

In future work, we plan on addressing these issues. Many of the observational uncertainties in the 
model are at high redshift. Modern 'semi-analytic' models {e.g., Cattaneo et al. 1999; Haehnelt & Kauff- 
mann 2000) go a long way toward elucidating these issues within the context of particular scenarios for 
structure formation. Observationally, the quasar luminosity function (Boyle et al. 2000; Fan ct al. 2001), 
OH Megamasers (Darling & Giovanelli 2002) and "X-type" radio sources (Merritt & Ekers 2002) are other 
tools for understanding the high-redshift MBH population. For this work, a detailed high-resolution search 
for more nearby MBH binaries themselves, perhaps in the form of double- nucleus QSOs (but see Kochanek 
et al. 1999) and Ultra-Luminous Infrared Galaxies (Murphy ct al. 2001), will be a crucial tool. Analytical 
and N-body studies of the dynamics and energetics of MBH binaries within galaxies will provide population 
information on MBH-MBH mass ratios and orbit sizes crucial to future work. These will allow us to also 
constrain the dynamics of the MBH binary population: if and when do they stall prior to the BW regime? 

We can, of course, use our current model to at least see the effect of these uncertainties. We can find 

a lower limit to the effect of MBH mass-function evolution by using a suitable initial mass function for the 
whole redshift range. If the average MBH mass has increased by an order of magnitude (say), then this 
would lower the characteristic strain (eq. [32]) by no more than 10^/*^ w 7. 

The current experimental limit on the low- frequency gravitational wave background has improved over 
the pioneering work of Kaspi et al. (1994), and will be accurately stated by Lommen & Backer (in prepa- 
ration). The new limit is lower at lower frequencies while the expected spcctnim is rising. This new result 
will constrain parameters of our model spectrum on the high side of the current estimate. Ongoing work 
with an array of millisecond pulsars has the prospect of significant improvement of detection capability in 
the coming decade. Techniques for optimal fitting of Pulsar Timing Array data need to be further developed 
to meet the demands of the new measurements. 

We thank Jon Arons, Ray Carlberg, Ron Hellings, Yuri Levin, Andrea Lommen, John Maggorian, Sterl 
Phinney, Philip Stark and members of the Center for Particle Astrophysics for helpful conversations. AHJ 
acknowledges support from NSF KDI grant 9872979 and NASA LTSA grant NAG5-6552, and PPARC in 
the UK. DCB acknowledges support from NSF grant AST-9731106 which partially supported the Pulsar 
Timing Array experiment. 



In this appendix, we calculate the probability distribution of the quantity h^{f) as given by equation (28). 
For generality, consider some quantity 



where N{z) dz gives the probability of some event happening in (z, z + dz) (i.e., a Poisson rate), and g{z) 
is any function. We want to know P{y), the distribution function of y. We start with a few important facts 
from probability theory. 



A. Appendix 




(Al) 



-31 - 



First, the characteristic function (or moment generating function) of a distribution is defined as the 
Fourier transform of the distribution, i.e., 

^{t) = {e'*y) =Jdy e'^ypiy) . (A2) 

This has the property that the moments of the distribution are then given by 

(A3) 



and the "cumulants" or "connected moments" or "semi- invariants" (k^., the mean, variance, skewness, kur- 
tosis, . . . , for k = 1, 2, 3, 4, . . .) are given by the same expression, substituting hup for ip. 

Next, we use the fact (trivial from Fourier theory) that the characteristic function of the sum of two 
variables is the product of the individual functions, and that the characteristic function of ay, where a is a 
scalar and y is the random variable is ip{at), where ip{t) is the characteristic function for y. 

Finally, we need the characteristic function for a Poisson distribution, P{n) = r"e~''/n! where r is the 
rate. This is just ip{t) = exp [r(e** — 1)] . 

Putting all of these facts together gives the characteristic function for y, which is just the sum (integral) 
of Poisson variables with rates N{z) dz times scalars g{z). Thus 

In ip{t) = J dz N{z) [e"s(^) _ ij . (A4) 

The Taylor series of In (p around t = then gives the cumulants of the distribution. They are just 

Kk = j dz [9{ztN{z) (A5) 

(Note that the mean, fc = 1, is what would be trivially expected.) 

In our case, we have N{z) dz 7V(z, /, Mi, M2) dz dMi dM2 and g{z) h'^^^{z, f, Mi, M2) for 
y = ^c(/) df/f, giving equation (36) above. 
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